home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 1 Issue 2 / PDCD-1 - Issue 02.iso / _utilities / utilities / 001 / meschach / !Meschach / c / vecop < prev    next >
Text File  |  1994-03-08  |  13KB  |  606 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /* vecop.c 1.3 8/18/87 */
  28.  
  29. #include    <stdio.h>
  30. #include    "matrix.h"
  31.  
  32. static    char    rcsid[] = "$Id: vecop.c,v 1.4 1994/03/08 05:50:39 des Exp $";
  33.  
  34.  
  35. /* _in_prod -- inner product of two vectors from i0 downwards */
  36. double    _in_prod(a,b,i0)
  37. VEC    *a,*b;
  38. u_int    i0;
  39. {
  40.     u_int    limit;
  41.     /* Real    *a_v, *b_v; */
  42.     /* register Real    sum; */
  43.  
  44.     if ( a==(VEC *)NULL || b==(VEC *)NULL )
  45.         error(E_NULL,"_in_prod");
  46.     limit = min(a->dim,b->dim);
  47.     if ( i0 > limit )
  48.         error(E_BOUNDS,"_in_prod");
  49.  
  50.     return __ip__(&(a->ve[i0]),&(b->ve[i0]),(int)(limit-i0));
  51.     /*****************************************
  52.     a_v = &(a->ve[i0]);        b_v = &(b->ve[i0]);
  53.     for ( i=i0; i<limit; i++ )
  54.         sum += a_v[i]*b_v[i];
  55.         sum += (*a_v++)*(*b_v++);
  56.  
  57.     return (double)sum;
  58.     ******************************************/
  59. }
  60.  
  61. /* sv_mlt -- scalar-vector multiply -- may be in-situ */
  62. VEC    *sv_mlt(scalar,vector,out)
  63. double    scalar;
  64. VEC    *vector,*out;
  65. {
  66.     /* u_int    dim, i; */
  67.     /* Real    *out_ve, *vec_ve; */
  68.  
  69.     if ( vector==(VEC *)NULL )
  70.         error(E_NULL,"sv_mlt");
  71.     if ( out==(VEC *)NULL || out->dim != vector->dim )
  72.         out = v_resize(out,vector->dim);
  73.     if ( scalar == 0.0 )
  74.         return v_zero(out);
  75.     if ( scalar == 1.0 )
  76.         return v_copy(vector,out);
  77.  
  78.     __smlt__(vector->ve,(double)scalar,out->ve,(int)(vector->dim));
  79.     /**************************************************
  80.     dim = vector->dim;
  81.     out_ve = out->ve;    vec_ve = vector->ve;
  82.     for ( i=0; i<dim; i++ )
  83.         out->ve[i] = scalar*vector->ve[i];
  84.         (*out_ve++) = scalar*(*vec_ve++);
  85.     **************************************************/
  86.     return (out);
  87. }
  88.  
  89. /* v_add -- vector addition -- may be in-situ */
  90. VEC    *v_add(vec1,vec2,out)
  91. VEC    *vec1,*vec2,*out;
  92. {
  93.     u_int    dim;
  94.     /* Real    *out_ve, *vec1_ve, *vec2_ve; */
  95.  
  96.     if ( vec1==(VEC *)NULL || vec2==(VEC *)NULL )
  97.         error(E_NULL,"v_add");
  98.     if ( vec1->dim != vec2->dim )
  99.         error(E_SIZES,"v_add");
  100.     if ( out==(VEC *)NULL || out->dim != vec1->dim )
  101.         out = v_resize(out,vec1->dim);
  102.     dim = vec1->dim;
  103.     __add__(vec1->ve,vec2->ve,out->ve,(int)dim);
  104.     /************************************************************
  105.     out_ve = out->ve;    vec1_ve = vec1->ve;    vec2_ve = vec2->ve;
  106.     for ( i=0; i<dim; i++ )
  107.         out->ve[i] = vec1->ve[i]+vec2->ve[i];
  108.         (*out_ve++) = (*vec1_ve++) + (*vec2_ve++);
  109.     ************************************************************/
  110.  
  111.     return (out);
  112. }
  113.  
  114. /* v_mltadd -- scalar/vector multiplication and addition
  115.         -- out = v1 + scale.v2        */
  116. VEC    *v_mltadd(v1,v2,scale,out)
  117. VEC    *v1,*v2,*out;
  118. double    scale;
  119. {
  120.     /* register u_int    dim, i; */
  121.     /* Real    *out_ve, *v1_ve, *v2_ve; */
  122.  
  123.     if ( v1==(VEC *)NULL || v2==(VEC *)NULL )
  124.         error(E_NULL,"v_mltadd");
  125.     if ( v1->dim != v2->dim )
  126.         error(E_SIZES,"v_mltadd");
  127.     if ( scale == 0.0 )
  128.         return v_copy(v1,out);
  129.     if ( scale == 1.0 )
  130.         return v_add(v1,v2,out);
  131.  
  132.     if ( v2 != out )
  133.     {
  134.         tracecatch(out = v_copy(v1,out),"v_mltadd");
  135.  
  136.         /* dim = v1->dim; */
  137.         __mltadd__(out->ve,v2->ve,scale,(int)(v1->dim));
  138.     }
  139.     else
  140.     {
  141.         tracecatch(out = sv_mlt(scale,v2,out),"v_mltadd");
  142.         out = v_add(v1,out,out);
  143.     }
  144.     /************************************************************
  145.     out_ve = out->ve;    v1_ve = v1->ve;        v2_ve = v2->ve;
  146.     for ( i=0; i < dim ; i++ )
  147.         out->ve[i] = v1->ve[i] + scale*v2->ve[i];
  148.         (*out_ve++) = (*v1_ve++) + scale*(*v2_ve++);
  149.     ************************************************************/
  150.  
  151.     return (out);
  152. }
  153.  
  154. /* v_sub -- vector subtraction -- may be in-situ */
  155. VEC    *v_sub(vec1,vec2,out)
  156. VEC    *vec1,*vec2,*out;
  157. {
  158.     /* u_int    i, dim; */
  159.     /* Real    *out_ve, *vec1_ve, *vec2_ve; */
  160.  
  161.     if ( vec1==(VEC *)NULL || vec2==(VEC *)NULL )
  162.         error(E_NULL,"v_sub");
  163.     if ( vec1->dim != vec2->dim )
  164.         error(E_SIZES,"v_sub");
  165.     if ( out==(VEC *)NULL || out->dim != vec1->dim )
  166.         out = v_resize(out,vec1->dim);
  167.  
  168.     __sub__(vec1->ve,vec2->ve,out->ve,(int)(vec1->dim));
  169.     /************************************************************
  170.     dim = vec1->dim;
  171.     out_ve = out->ve;    vec1_ve = vec1->ve;    vec2_ve = vec2->ve;
  172.     for ( i=0; i<dim; i++ )
  173.         out->ve[i] = vec1->ve[i]-vec2->ve[i];
  174.         (*out_ve++) = (*vec1_ve++) - (*vec2_ve++);
  175.     ************************************************************/
  176.  
  177.     return (out);
  178. }
  179.  
  180. /* v_map -- maps function f over components of x: out[i] = f(x[i])
  181.     -- _v_map sets out[i] = f(params,x[i]) */
  182. VEC    *v_map(f,x,out)
  183. #ifdef PROTOTYPES_IN_STRUCT
  184. double    (*f)(double);
  185. #else
  186. double    (*f)();
  187. #endif
  188. VEC    *x, *out;
  189. {
  190.     Real    *x_ve, *out_ve;
  191.     int    i, dim;
  192.  
  193.     if ( ! x || ! f )
  194.         error(E_NULL,"v_map");
  195.     if ( ! out || out->dim != x->dim )
  196.         out = v_resize(out,x->dim);
  197.  
  198.     dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  199.     for ( i = 0; i < dim; i++ )
  200.         *out_ve++ = (*f)(*x_ve++);
  201.  
  202.     return out;
  203. }
  204.  
  205. VEC    *_v_map(f,params,x,out)
  206. #ifdef PROTOTYPES_IN_STRUCT
  207. double    (*f)(void *,double);
  208. #else
  209. double    (*f)();
  210. #endif
  211. VEC    *x, *out;
  212. void    *params;
  213. {
  214.     Real    *x_ve, *out_ve;
  215.     int    i, dim;
  216.  
  217.     if ( ! x || ! f )
  218.         error(E_NULL,"_v_map");
  219.     if ( ! out || out->dim != x->dim )
  220.         out = v_resize(out,x->dim);
  221.  
  222.     dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  223.     for ( i = 0; i < dim; i++ )
  224.         *out_ve++ = (*f)(params,*x_ve++);
  225.  
  226.     return out;
  227. }
  228.  
  229. /* v_lincomb -- returns sum_i a[i].v[i], a[i] real, v[i] vectors */
  230. VEC    *v_lincomb(n,v,a,out)
  231. int    n;    /* number of a's and v's */
  232. Real    a[];
  233. VEC    *v[], *out;
  234. {
  235.     int    i;
  236.  
  237.     if ( ! a || ! v )
  238.         error(E_NULL,"v_lincomb");
  239.     if ( n <= 0 )
  240.         return VNULL;
  241.  
  242.     for ( i = 1; i < n; i++ )
  243.         if ( out == v[i] )
  244.             error(E_INSITU,"v_lincomb");
  245.  
  246.     out = sv_mlt(a[0],v[0],out);
  247.     for ( i = 1; i < n; i++ )
  248.     {
  249.         if ( ! v[i] )
  250.             error(E_NULL,"v_lincomb");
  251.         if ( v[i]->dim != out->dim )
  252.             error(E_SIZES,"v_lincomb");
  253.         out = v_mltadd(out,v[i],a[i],out);
  254.     }
  255.  
  256.     return out;
  257. }
  258.  
  259.  
  260.  
  261. #ifdef ANSI_C
  262.  
  263. /* v_linlist -- linear combinations taken from a list of arguments;
  264.    calling:
  265.       v_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  266.    where vi are vectors (VEC *) and ai are numbers (double)
  267. */
  268. VEC  *v_linlist(VEC *out,VEC *v1,double a1,...)
  269. {
  270.    va_list ap;
  271.    VEC *par;
  272.    double a_par;
  273.  
  274.    if ( ! v1 )
  275.      return VNULL;
  276.    
  277.    va_start(ap, a1);
  278.    out = sv_mlt(a1,v1,out);
  279.    
  280.    while (par = va_arg(ap,VEC *)) {   /* NULL ends the list*/
  281.       a_par = va_arg(ap,double);
  282.       if (a_par == 0.0) continue;
  283.       if ( out == par )        
  284.     error(E_INSITU,"v_linlist");
  285.       if ( out->dim != par->dim )    
  286.     error(E_SIZES,"v_linlist");
  287.  
  288.       if (a_par == 1.0)
  289.     out = v_add(out,par,out);
  290.       else if (a_par == -1.0)
  291.     out = v_sub(out,par,out);
  292.       else
  293.     out = v_mltadd(out,par,a_par,out); 
  294.    } 
  295.    
  296.    va_end(ap);
  297.    return out;
  298. }
  299.  
  300. #elif VARARGS
  301.  
  302.  
  303. /* v_linlist -- linear combinations taken from a list of arguments;
  304.    calling:
  305.       v_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  306.    where vi are vectors (VEC *) and ai are numbers (double)
  307. */
  308. VEC  *v_linlist(va_alist) va_dcl
  309. {
  310.    va_list ap;
  311.    VEC *par, *out;
  312.    double a_par;
  313.  
  314.    va_start(ap);
  315.    out = va_arg(ap,VEC *);
  316.    par = va_arg(ap,VEC *);
  317.    if ( ! par ) {
  318.       va_end(ap);
  319.       return VNULL;
  320.    }
  321.    
  322.    a_par = va_arg(ap,double);
  323.    out = sv_mlt(a_par,par,out);
  324.    
  325.    while (par = va_arg(ap,VEC *)) {   /* NULL ends the list*/
  326.       a_par = va_arg(ap,double);
  327.       if (a_par == 0.0) continue;
  328.       if ( out == par )        
  329.     error(E_INSITU,"v_linlist");
  330.       if ( out->dim != par->dim )    
  331.     error(E_SIZES,"v_linlist");
  332.  
  333.       if (a_par == 1.0)
  334.     out = v_add(out,par,out);
  335.       else if (a_par == -1.0)
  336.     out = v_sub(out,par,out);
  337.       else
  338.     out = v_mltadd(out,par,a_par,out); 
  339.    } 
  340.    
  341.    va_end(ap);
  342.    return out;
  343. }
  344.  
  345. #endif
  346.   
  347.  
  348.  
  349.  
  350.  
  351. /* v_star -- computes componentwise (Hadamard) product of x1 and x2
  352.     -- result out is returned */
  353. VEC    *v_star(x1, x2, out)
  354. VEC    *x1, *x2, *out;
  355. {
  356.     int        i;
  357.  
  358.     if ( ! x1 || ! x2 )
  359.     error(E_NULL,"v_star");
  360.     if ( x1->dim != x2->dim )
  361.     error(E_SIZES,"v_star");
  362.     out = v_resize(out,x1->dim);
  363.  
  364.     for ( i = 0; i < x1->dim; i++ )
  365.     out->ve[i] = x1->ve[i] * x2->ve[i];
  366.  
  367.     return out;
  368. }
  369.  
  370. /* v_slash -- computes componentwise ratio of x2 and x1
  371.     -- out[i] = x2[i] / x1[i]
  372.     -- if x1[i] == 0 for some i, then raise E_SING error
  373.     -- result out is returned */
  374. VEC    *v_slash(x1, x2, out)
  375. VEC    *x1, *x2, *out;
  376. {
  377.     int        i;
  378.     Real    tmp;
  379.  
  380.     if ( ! x1 || ! x2 )
  381.     error(E_NULL,"v_slash");
  382.     if ( x1->dim != x2->dim )
  383.     error(E_SIZES,"v_slash");
  384.     out = v_resize(out,x1->dim);
  385.  
  386.     for ( i = 0; i < x1->dim; i++ )
  387.     {
  388.     tmp = x1->ve[i];
  389.     if ( tmp == 0.0 )
  390.         error(E_SING,"v_slash");
  391.     out->ve[i] = x2->ve[i] / tmp;
  392.     }
  393.  
  394.     return out;
  395. }
  396.  
  397. /* v_min -- computes minimum component of x, which is returned
  398.     -- also sets min_idx to the index of this minimum */
  399. double    v_min(x, min_idx)
  400. VEC    *x;
  401. int    *min_idx;
  402. {
  403.     int        i, i_min;
  404.     Real    min_val, tmp;
  405.  
  406.     if ( ! x )
  407.     error(E_NULL,"v_min");
  408.     if ( x->dim <= 0 )
  409.     error(E_SIZES,"v_min");
  410.     i_min = 0;
  411.     min_val = x->ve[0];
  412.     for ( i = 1; i < x->dim; i++ )
  413.     {
  414.     tmp = x->ve[i];
  415.     if ( tmp < min_val )
  416.     {
  417.         min_val = tmp;
  418.         i_min = i;
  419.     }
  420.     }
  421.  
  422.     if ( min_idx != NULL )
  423.     *min_idx = i_min;
  424.     return min_val;
  425. }
  426.  
  427. /* v_max -- computes maximum component of x, which is returned
  428.     -- also sets max_idx to the index of this maximum */
  429. double    v_max(x, max_idx)
  430. VEC    *x;
  431. int    *max_idx;
  432. {
  433.     int        i, i_max;
  434.     Real    max_val, tmp;
  435.  
  436.     if ( ! x )
  437.     error(E_NULL,"v_max");
  438.     if ( x->dim <= 0 )
  439.     error(E_SIZES,"v_max");
  440.     i_max = 0;
  441.     max_val = x->ve[0];
  442.     for ( i = 1; i < x->dim; i++ )
  443.     {
  444.     tmp = x->ve[i];
  445.     if ( tmp > max_val )
  446.     {
  447.         max_val = tmp;
  448.         i_max = i;
  449.     }
  450.     }
  451.  
  452.     if ( max_idx != NULL )
  453.     *max_idx = i_max;
  454.     return max_val;
  455. }
  456.  
  457. #define    MAX_STACK    60
  458.  
  459.  
  460. /* v_sort -- sorts vector x, and generates permutation that gives the order
  461.     of the components; x = [1.3, 3.7, 0.5] -> [0.5, 1.3, 3.7] and
  462.     the permutation is order = [2, 0, 1].
  463.     -- if order is NULL on entry then it is ignored
  464.     -- the sorted vector x is returned */
  465. VEC    *v_sort(x, order)
  466. VEC    *x;
  467. PERM    *order;
  468. {
  469.     Real    *x_ve, tmp, v;
  470.     /* int        *order_pe; */
  471.     int        dim, i, j, l, r, tmp_i;
  472.     int        stack[MAX_STACK], sp;
  473.  
  474.     if ( ! x )
  475.     error(E_NULL,"v_sort");
  476.     if ( order != PNULL && order->size != x->dim )
  477.     order = px_resize(order, x->dim);
  478.  
  479.     x_ve = x->ve;
  480.     dim = x->dim;
  481.     if ( order != PNULL )
  482.     px_ident(order);
  483.  
  484.     if ( dim <= 1 )
  485.     return x;
  486.  
  487.     /* using quicksort algorithm in Sedgewick,
  488.        "Algorithms in C", Ch. 9, pp. 118--122 (1990) */
  489.     sp = 0;
  490.     l = 0;    r = dim-1;    v = x_ve[0];
  491.     for ( ; ; )
  492.     {
  493.     while ( r > l )
  494.     {
  495.         /* "i = partition(x_ve,l,r);" */
  496.         v = x_ve[r];
  497.         i = l-1;
  498.         j = r;
  499.         for ( ; ; )
  500.         {
  501.         while ( x_ve[++i] < v )
  502.             ;
  503.         while ( x_ve[--j] > v )
  504.             ;
  505.         if ( i >= j )    break;
  506.         
  507.         tmp = x_ve[i];
  508.         x_ve[i] = x_ve[j];
  509.         x_ve[j] = tmp;
  510.         if ( order != PNULL )
  511.         {
  512.             tmp_i = order->pe[i];
  513.             order->pe[i] = order->pe[j];
  514.             order->pe[j] = tmp_i;
  515.         }
  516.         }
  517.         tmp = x_ve[i];
  518.         x_ve[i] = x_ve[r];
  519.         x_ve[r] = tmp;
  520.         if ( order != PNULL )
  521.         {
  522.         tmp_i = order->pe[i];
  523.         order->pe[i] = order->pe[r];
  524.         order->pe[r] = tmp_i;
  525.         }
  526.  
  527.         if ( i-l > r-i )
  528.         {   stack[sp++] = l;   stack[sp++] = i-1;   l = i+1;   }
  529.         else
  530.         {   stack[sp++] = i+1;   stack[sp++] = r;   r = i-1;   }
  531.     }
  532.  
  533.     /* recursion elimination */
  534.     if ( sp == 0 )
  535.         break;
  536.     r = stack[--sp];
  537.     l = stack[--sp];
  538.     }
  539.  
  540.     return x;
  541. }
  542.  
  543. /* v_sum -- returns sum of entries of a vector */
  544. double    v_sum(x)
  545. VEC    *x;
  546. {
  547.     int        i;
  548.     Real    sum;
  549.  
  550.     if ( ! x )
  551.     error(E_NULL,"v_sum");
  552.  
  553.     sum = 0.0;
  554.     for ( i = 0; i < x->dim; i++ )
  555.     sum += x->ve[i];
  556.  
  557.     return sum;
  558. }
  559.  
  560. /* v_conv -- computes convolution product of two vectors */
  561. VEC    *v_conv(x1, x2, out)
  562. VEC    *x1, *x2, *out;
  563. {
  564.     int        i;
  565.  
  566.     if ( ! x1 || ! x2 )
  567.     error(E_NULL,"v_conv");
  568.     if ( x1 == out || x2 == out )
  569.     error(E_INSITU,"v_conv");
  570.     if ( x1->dim == 0 || x2->dim == 0 )
  571.     return out = v_resize(out,0);
  572.  
  573.     out = v_resize(out,x1->dim + x2->dim - 1);
  574.     v_zero(out);
  575.     for ( i = 0; i < x1->dim; i++ )
  576.     __mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim);
  577.  
  578.     return out;
  579. }
  580.  
  581. /* v_pconv -- computes a periodic convolution product
  582.     -- the period is the dimension of x2 */
  583. VEC    *v_pconv(x1, x2, out)
  584. VEC    *x1, *x2, *out;
  585. {
  586.     int        i;
  587.  
  588.     if ( ! x1 || ! x2 )
  589.     error(E_NULL,"v_pconv");
  590.     if ( x1 == out || x2 == out )
  591.     error(E_INSITU,"v_pconv");
  592.     out = v_resize(out,x2->dim);
  593.     if ( x2->dim == 0 )
  594.     return out;
  595.  
  596.     v_zero(out);
  597.     for ( i = 0; i < x1->dim; i++ )
  598.     {
  599.     __mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim - i);
  600.     if ( i > 0 )
  601.         __mltadd__(out->ve,&(x2->ve[x2->dim - i]),x1->ve[i],i);
  602.     }
  603.  
  604.     return out;
  605. }
  606.